
#include "timeDiscr.h"
#include <iostream>
#include <fstream>
using namespace std;

TimeDiscr::TimeDiscr(const Config &cfg, const Geometry &geo, Fluid &flow,
                     SpaceDiscr &spaceDiscrete, BoundConds &bound)
    : cv{flow.GetCV()},
      si{geo.si}, sj{geo.sj}, vol{geo.vol},
      dvE{flow.GetdvEuler()}, dvNS{flow.GetdvNS()}, coord{geo.coord},
      rhs{spaceDiscrete.GetRhs()},
      timeStepType{cfg.GetTimeStepType()}, numStage{cfg.GetNumStage()},
      cflNum{cfg.GetCflNum()}, ark{cfg.GetStageCoef()},
      timeDiscType{cfg.GetTimeDiscType()},
      equType{cfg.GetEquationType()},
      fluxType{cfg.GetConvectFluxType()},
      spaceDiscrt{spaceDiscrete},
      flows{flow}, boundary{bound}, imax{geo.imax},
      jmax{geo.jmax}, il{imax - 1}, jl{jmax - 1}, i2{il - 1}, j2{jl - 1}

{
  int imaxp1{imax + 1}, jmaxp1{jmax + 1};
  cvOld.ReConstruct(imaxp1, jmaxp1);
  tstep.ReConstruct(imaxp1, jmaxp1);
}

void TimeDiscr::SolveOneStepEuler(void)
{
  cvOld = cv;
  TimeStepEuler();
  real_t fac, adtv;
  for (int stage = 0; stage < numStage; ++stage)
  {
    spaceDiscrt.SetDissipationZero();

    spaceDiscrt.DissipRoe1();
    spaceDiscrt.FluxRoe1();

    fac = ark[stage] * cflNum;
    for (int j = 2; j <= j2; ++j)
    {
      for (int i = 2; i <= i2; ++i)
      {
        adtv = fac * tstep(i, j) / vol(i, j);
        rhs(i, j).rho *= adtv;
        rhs(i, j).rhou *= adtv;
        rhs(i, j).rhov *= adtv;
        rhs(i, j).rhoE *= adtv;
      }
    }
    // implicit residual smoothing   all smoothing
    // 有待添加

    // update - new solution, new dependent variables
    for (int j = 2; j <= j2; ++j)
    {
      for (int i = 2; i <= i2; ++i)
      {
        cv(i, j).rho = cvOld(i, j).rho - rhs(i, j).rho;
        cv(i, j).rhou = cvOld(i, j).rhou - rhs(i, j).rhou;
        cv(i, j).rhov = cvOld(i, j).rhov - rhs(i, j).rhov;
        cv(i, j).rhoE = cvOld(i, j).rhoE - rhs(i, j).rhoE;
      }
    }

    flows.DependentVarsAllEuler();
    boundary.BoundaryConditionsEuler();
  } // end irk
}

//===========================================
void TimeDiscr::SolveOneStepNS(void)
{
  cvOld = cv;
  TimeStepNS();
  real_t fac, adtv;
  for (int stage = 0; stage < numStage; ++stage)
  {
    spaceDiscrt.SetDissipationZero();

    spaceDiscrt.GradientInitial();

    // for (int j = 1; j <= jmax; ++j)
    // {
    //   for (int i = 1; i <= imax; ++i)
    //   {
    //     spaceDiscrt.u(i, j) = (i + j) * 1.0;
    //     spaceDiscrt.v(i, j) = (i + j) * 1.0;
    //     dvNS(i, j).T = (i + j) * 1.0;
    //     dvNS(i, j).ka = (i + j) * 1.0;
    //     dvNS(i, j).mu = (i + j) * 1.0;
    //   }
    // }

    spaceDiscrt.GradientFaceI(); // 没有问题

    // std::ofstream gradfif{"checkGradfi"}, gradfjf{"checkGradfj"};
    // printCheck(gradfif, spaceDiscrt.gradFi, 0, imax, 0, jmax);

    spaceDiscrt.GradientFaceJ(); // 没有问题

    // printCheck(gradfjf, spaceDiscrt.gradFj, 0, imax, 0, jmax);
    // gradfif.flush();
    // gradfjf.flush();
    // spaceDiscrt.SetDissipationZero();

    spaceDiscrt.FluxViscous(); // 没问题

    // std::ofstream fluxviscousf{"checkfluxviscous"};
    // printCheck(fluxviscousf, spaceDiscrt.diss, 0, imax, 0, jmax);
    // fluxviscousf.flush();
    // exit(EXIT_FAILURE);

    spaceDiscrt.DissipRoe1();
    spaceDiscrt.FluxRoe1();

    fac = ark[stage] * cflNum;
    for (int j = 2; j <= j2; ++j)
    {
      for (int i = 2; i <= i2; ++i)
      {
        adtv = fac * tstep(i, j) / vol(i, j);
        rhs(i, j).rho *= adtv;
        rhs(i, j).rhou *= adtv;
        rhs(i, j).rhov *= adtv;
        rhs(i, j).rhoE *= adtv;
      }
    }
    // implicit residual smoothing   all smoothing
    // 有待添加

    // update - new solution, new dependent variables
    for (int j = 2; j <= j2; ++j)
    {
      for (int i = 2; i <= i2; ++i)
      {
        cv(i, j).rho = cvOld(i, j).rho - rhs(i, j).rho;
        cv(i, j).rhou = cvOld(i, j).rhou - rhs(i, j).rhou;
        cv(i, j).rhov = cvOld(i, j).rhov - rhs(i, j).rhov;
        cv(i, j).rhoE = cvOld(i, j).rhoE - rhs(i, j).rhoE;
      }
    }

    flows.DependVarsAllNS();
    boundary.BoundaryConditionsNS();
  } // end irk
}

void TimeDiscr::DensityChange2Norm(real_t &drho2Norm, real_t &drhoMax,
                                   Coord_t &position)
{
  real_t dr;
  drho2Norm = drhoMax = 0.0;
  for (int j = 2; j <= j2; ++j)
  {
    for (int i = 2; i <= i2; ++i)
    {
      dr = cv(i, j).rho - cvOld(i, j).rho;
      // cout << dr << endl;
      drho2Norm += dr * dr;
      if (drhoMax < fabs(dr))
      {
        drhoMax = fabs(dr);
        position.x = 0.25 * (coord(i, j).x + coord(i + 1, j).x +
                             coord(i, j + 1).x + coord(i + 1, j + 1).x);
        position.y = 0.25 * (coord(i, j).y + coord(i + 1, j).y +
                             coord(i, j + 1).y + coord(i + 1, j + 1).y);
      }
    }
  }
  drho2Norm = sqrt(drho2Norm);
}

void TimeDiscr::TimeStepEuler(void)
{
  real_t rrho, u, v, sx, sy, vc, cs, sri, srj;

  for (int j = 2; j <= j2; ++j)
  {
    for (int i = 2; i <= i2; ++i)
    {
      rrho = 1.0 / cv(i, j).rho;
      u = cv(i, j).rhou * rrho;
      v = cv(i, j).rhov * rrho;

      sx = 0.5 * (si(i, j).sx + si(i + 1, j).sx);
      sy = 0.5 * (si(i, j).sy + si(i + 1, j).sy);
      vc = sx * u + sy * v;
      cs = dvE(i, j).c * sqrt(sx * sx + sy * sy);
      sri = fabs(vc) + cs;

      sx = 0.5 * (sj(i, j).sx + sj(i, j + 1).sx);
      sy = 0.5 * (sj(i, j).sy + sj(i, j + 1).sy);
      vc = sx * u + sy * v;
      cs = dvE(i, j).c * sqrt(sx * sx + sy * sy);
      srj = fabs(vc) + cs;

      tstep(i, j) = vol(i, j) / (sri + srj);
    }
  }

  if (timeStepType == TimeStep_t::Global)
  {
    real_t tmin{1.0e10};
    for (int j = 2; j <= j2; ++j)
      for (int i = 2; i <= i2; ++i)
        tmin = fmin(tstep(i, j), tmin);

    tstep = tmin;
  }
}

void TimeDiscr::TimeStepNS(void)
{
  real_t rrho, u, v, sx, sy, vc, cs, sri,
      srj, cfac, coe, ds2, srvi, srvj;

  // 后面有可能添加通量格式，这个保险  会有提示
  switch (fluxType)
  {
  case ConvectFlux_t::Roe1:
    cfac = 2.0;
    break;

  case ConvectFlux_t::Roe2:
    cfac = 1.0;
    break;
  }

  for (int j = 2; j <= j2; ++j)
  {
    for (int i = 2; i <= i2; ++i)
    {
      rrho = 1.0 / cv(i, j).rho;
      u = cv(i, j).rhou * rrho;
      v = cv(i, j).rhov * rrho;
      coe = gama * rrho * dvNS(i, j).mu / Prl / vol(i, j);

      sx = 0.5 * (si(i, j).sx + si(i + 1, j).sx);
      sy = 0.5 * (si(i, j).sy + si(i + 1, j).sy);
      ds2 = sx * sx + sy * sy;
      vc = sx * u + sy * v;
      cs = dvNS(i, j).c * sqrt(ds2);
      sri = fabs(vc) + cs;
      srvi = coe * ds2;

      sx = 0.5 * (sj(i, j).sx + sj(i, j + 1).sx);
      sy = 0.5 * (sj(i, j).sy + sj(i, j + 1).sy);
      ds2 = sx * sx + sy * sy;
      vc = sx * u + sy * v;
      cs = dvNS(i, j).c * sqrt(ds2);
      srj = fabs(vc) + cs;
      srvj = coe * ds2;

      tstep(i, j) = vol(i, j) / (sri + srj + cfac * (srvi + srvj));
    }
  }

  if (timeStepType == TimeStep_t::Global)
  {
    real_t tmin{1.0e10};
    for (int j = 2; j <= j2; ++j)
      for (int i = 2; i <= i2; ++i)
        tmin = fmin(tstep(i, j), tmin);

    tstep = tmin;
  }
}
